시계열분석 8장

library(glue)
list.files("../adv_reg/functions")

## [1] "anova_reg.R"       "biased_reg.R"      "chapter_11"       
## [4] "mult_reg.R"        "reg_diagnostics.R" "var_selection.R"

dir = getwd() ; mother_dir = dirname(dir)
dataset_dir = glue(dir, "/datasets", sep = "")

reg_function_dir = glue(mother_dir, "/adv_reg/functions")
source(glue(reg_function_dir, "/mult_reg.R", sep = ""), echo = F)
source(glue(reg_function_dir, "/reg_diagnostics.R", sep = ""), echo = F)

tsa_function_dir = glue(mother_dir, "/tsa/functions")
source(glue(tsa_function_dir, "/exp_smt.R", sep = ""), echo = F)
source(glue(tsa_function_dir, "/acf_pacf.R", sep = ""), echo = F)

8.5 모형의 적합 예제

(예 8-6)

## Example 8.6 : 가스로 자료 분석
# install.packages("portes") ; install.packages("lmtest")
library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(lmtest) # library for function coeftest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

gas = scan(glue(dataset_dir, "/gas.txt"), what = list(0,0))
names(gas) = c("rate", "co2") # 변수명 생성
gasrate = data.frame(gas)
time = 1:nrow(gasrate)
rate = ts(gas$rate)
co2 = ts(gas$co2)

ts.plot(rate, ylab = "gas rate", main = "그림 8-1 가스 공급비율")

acf2(rate, max.lag = 24, main = "그림 8-1 가스 공급비율")

##      [,1]  [,2] [,3] [,4] [,5]  [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.95  0.83 0.68 0.53 0.41  0.32 0.26 0.23 0.21  0.21  0.20  0.19  0.17
## PACF 0.95 -0.79 0.34 0.12 0.06 -0.11 0.05 0.10 0.02 -0.07 -0.09  0.04  0.09
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.14  0.10  0.08  0.05  0.04  0.03  0.04  0.06  0.07  0.08  0.08
## PACF -0.14  0.05  0.03 -0.02  0.03  0.09 -0.04 -0.09  0.04  0.04  0.01

fit1 = arima(rate, order = c(3,0,0)) # 절편 있는 AR(3) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        1.969066   0.054385  36.2061 < 2.2e-16 ***
## ar2       -1.365143   0.098538 -13.8540 < 2.2e-16 ***
## ar3        0.339404   0.054328   6.2473 4.177e-10 ***
## intercept -0.060643   0.189800  -0.3195    0.7493    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit2 = arima(rate, order = c(3,0,0), include.mean = F) # 절편 없는 AR(3) 모형 적합
summary(fit2)

## 
## Call:
## arima(x = rate, order = c(3, 0, 0), include.mean = F)
## 
## Coefficients:
##          ar1      ar2     ar3
##       1.9696  -1.3659  0.3399
## s.e.  0.0544   0.0985  0.0543
## 
## sigma^2 estimated as 0.03531:  log likelihood = 72.52,  aic = -137.04
## 
## Training set error measures:
##                        ME      RMSE       MAE MPE MAPE      MASE        ACF1
## Training set -0.003418121 0.1879034 0.1305995 NaN  Inf 0.5134245 -0.03583041

ts.plot(fit2$residuals, main = "그림 8-3 잔차 시계열그림"); abline(h = 0)

acf2(fit2$residuals, max.lag = 24, main = "그림 8-4 잔차 SACF와 SPACF")

##       [,1] [,2] [,3]  [,4]  [,5] [,6] [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.04 0.07 0.06 -0.14 -0.01 0.06 0.01  0.00 -0.05  0.04  0.14 -0.08  0.10
## PACF -0.04 0.07 0.06 -0.15 -0.03 0.08 0.04 -0.03 -0.08  0.05  0.18 -0.08  0.03
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.04 -0.08  0.02  0.06 -0.05 -0.08  0.02  0.01  0.03  0.04  0.00
## PACF  0.06 -0.02 -0.03  0.06 -0.02 -0.11  0.02  0.05  0.03  0.01 -0.04

qqnorm(resid(fit2), main = "그림 8-5 잔차의 정규성검정")
qqline(resid(fit2), col = "red")

LjungBox(fit2, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df    p-value
##     6  10.21058  3 0.01685835
##    12  19.76879  9 0.01939412
##    18  27.70650 15 0.02348012
##    24  30.86553 21 0.07592519

sarima(rate, 3,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 0.073551 
## iter   2 value -0.381366
## iter   3 value -0.607123
## iter   4 value -1.135181
## iter   5 value -1.188992
## iter   6 value -1.311157
## iter   7 value -1.364387
## iter   8 value -1.458785
## iter   9 value -1.553747
## iter  10 value -1.643806
## iter  11 value -1.666390
## iter  12 value -1.667454
## iter  13 value -1.667491
## iter  14 value -1.667506
## iter  15 value -1.667509
## iter  16 value -1.667514
## iter  17 value -1.667518
## iter  18 value -1.667519
## iter  19 value -1.667519
## iter  20 value -1.667519
## iter  21 value -1.667519
## iter  22 value -1.667519
## iter  23 value -1.667519
## iter  23 value -1.667519
## iter  23 value -1.667519
## final  value -1.667519 
## converged
## initial  value -1.664075 
## iter   2 value -1.664100
## iter   3 value -1.664101
## iter   4 value -1.664102
## iter   5 value -1.664102
## iter   6 value -1.664103
## iter   7 value -1.664103
## iter   8 value -1.664103
## iter   9 value -1.664104
## iter  10 value -1.664104
## iter  10 value -1.664104
## final  value -1.664104 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1     1.9691 0.0544  36.2061  0.0000
## ar2    -1.3651 0.0985 -13.8540  0.0000
## ar3     0.3394 0.0543   6.2473  0.0000
## xmean  -0.0606 0.1898  -0.3195  0.7496
## 
## sigma^2 estimated as 0.03529587 on 292 degrees of freedom 
##  
## AIC = -0.4565466  AICc = -0.4560823  BIC = -0.3942095 
## 

sarima.for(rate, 6, 3, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 297 
## End = 302 
## Frequency = 1 
## [1] -0.26510614 -0.22955238 -0.18245672 -0.13931247 -0.10658382 -0.08505255
## 
## $se
## Time Series:
## Start = 297 
## End = 302 
## Frequency = 1 
## [1] 0.1878719 0.4149044 0.6283958 0.7956610 0.9103261 0.9807042

sarima.for(rate, 12, 3, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 297 
## End = 308 
## Frequency = 1 
##  [1] -0.26510614 -0.22955238 -0.18245672 -0.13931247 -0.10658382 -0.08505255
##  [7] -0.07269199 -0.06663826 -0.06428419 -0.06371784 -0.06376165 -0.06382208
## 
## $se
## Time Series:
## Start = 297 
## End = 308 
## Frequency = 1 
##  [1] 0.1878719 0.4149044 0.6283958 0.7956610 0.9103261 0.9807042 1.0199474
##  [8] 1.0401144 1.0498515 1.0543918 1.0565120 1.0575442

sarima.for(rate, 18, 3, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 297 
## End = 314 
## Frequency = 1 
##  [1] -0.26510614 -0.22955238 -0.18245672 -0.13931247 -0.10658382 -0.08505255
##  [7] -0.07269199 -0.06663826 -0.06428419 -0.06371784 -0.06376165 -0.06382208
## [13] -0.06368903 -0.06335944 -0.06291258 -0.06243747 -0.06200011 -0.06163585
## 
## $se
## Time Series:
## Start = 297 
## End = 314 
## Frequency = 1 
##  [1] 0.1878719 0.4149044 0.6283958 0.7956610 0.9103261 0.9807042 1.0199474
##  [8] 1.0401144 1.0498515 1.0543918 1.0565120 1.0575442 1.0580851 1.0583935
## [15] 1.0585820 1.0587021 1.0587793 1.0588285

sarima.for(rate, 24, 3, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 297 
## End = 320 
## Frequency = 1 
##  [1] -0.26510614 -0.22955238 -0.18245672 -0.13931247 -0.10658382 -0.08505255
##  [7] -0.07269199 -0.06663826 -0.06428419 -0.06371784 -0.06376165 -0.06382208
## [13] -0.06368903 -0.06335944 -0.06291258 -0.06243747 -0.06200011 -0.06163585
## [19] -0.06135439 -0.06114902 -0.06100522 -0.06090691 -0.06083993 -0.06079344
## 
## $se
## Time Series:
## Start = 297 
## End = 320 
## Frequency = 1 
##  [1] 0.1878719 0.4149044 0.6283958 0.7956610 0.9103261 0.9807042 1.0199474
##  [8] 1.0401144 1.0498515 1.0543918 1.0565120 1.0575442 1.0580851 1.0583935
## [15] 1.0585820 1.0587021 1.0587793 1.0588285 1.0588590 1.0588774 1.0588882
## [22] 1.0588944 1.0588979 1.0588998

(예 8-7)

## Example 8.7 : 과대적합
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest

z = scan(glue("{dataset_dir}/eg8_7.txt"))
z.ts = ts(z)
ts.plot(z.ts, ylab = "z", main = "그림 8-6 모의실험 자료")

acf2(z.ts, max.lag = 24, main = "그림 8-6 ACF & PACF")

##       [,1] [,2]  [,3] [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.68 0.50 -0.34 0.24 -0.15  0.07 -0.12  0.05 -0.10  0.10 -0.19  0.17
## PACF -0.68 0.08  0.04 0.00  0.02 -0.06 -0.18 -0.12 -0.13  0.01 -0.16 -0.07
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.11  0.14 -0.14  0.18 -0.16  0.12 -0.09  0.04 -0.11  0.08 -0.10  0.09
## PACF  0.07  0.10 -0.07  0.05 -0.03 -0.13 -0.04 -0.02 -0.17 -0.10 -0.05  0.01

LjungBox(z.ts, lag = seq(6,24,6))

##  lags statistic df      p-value
##     6  93.92521  6 0.000000e+00
##    12 105.36773 12 0.000000e+00
##    18 120.45361 18 0.000000e+00
##    24 126.81112 24 5.551115e-16

fit = arima(z.ts, order = c(1,0,0));
summary(fit)

## 
## Call:
## arima(x = z.ts, order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.6715    19.8312
## s.e.   0.0728     0.1776
## 
## sigma^2 estimated as 8.744:  log likelihood = -250.61,  aic = 507.22
## 
## Training set error measures:
##                       ME     RMSE      MAE       MPE    MAPE      MASE
## Training set 0.007670341 2.957007 2.337302 -2.080783 11.9899 0.4096807
##                    ACF1
## Training set 0.04344908

coeftest(fit)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1       -0.671544   0.072838  -9.2197 < 2.2e-16 ***
## intercept 19.831150   0.177618 111.6504 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ts.plot(fit$residuals, main = "그림 8-7 잔차", ylab = "residual") ; abline(h = 0)

acf2(fit$residuals, max.lag = 24, main = "그림 8-8 잔차 SACF와 SPACF")

##      [,1] [,2]  [,3] [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.04 0.07  0.00 0.03 -0.03 -0.14 -0.16 -0.12 -0.08 -0.08 -0.16  0.09  0.11
## PACF 0.04 0.07 -0.01 0.03 -0.03 -0.15 -0.15 -0.10 -0.06 -0.06 -0.16  0.09  0.09
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.07  0.02  0.11 -0.07  0.00 -0.04 -0.13 -0.15 -0.05 -0.06 -0.09
## PACF  0.01 -0.03  0.05 -0.16 -0.05 -0.02 -0.13 -0.15 -0.02 -0.02 -0.11

qqnorm(fit$residuals, main = "그림 8-9 잔차의 정규성검정")
qqline(fit$residuals, col = "red")

LjungBox(fit, lag = seq(6, 24, 6)) # 잔차의 포트맨토검정

##  lags statistic df   p-value
##     6  3.147291  5 0.6772898
##    12 13.135047 11 0.2845906
##    18 17.072699 17 0.4494514
##    24 24.393987 23 0.3822679

sarima.for(z.ts, 25, 1, 0, 0)

## $pred
## Time Series:
## Start = 101 
## End = 125 
## Frequency = 1 
##  [1] 20.82225 19.16558 20.27811 19.53100 20.03271 19.69579 19.92205 19.77011
##  [9] 19.87214 19.80362 19.84964 19.81874 19.83949 19.82555 19.83491 19.82863
## [17] 19.83285 19.83001 19.83191 19.83064 19.83149 19.83092 19.83131 19.83105
## [25] 19.83122
## 
## $se
## Time Series:
## Start = 101 
## End = 125 
## Frequency = 1 
##  [1] 2.957007 3.561901 3.803344 3.907350 3.953359 3.973933 3.983176 3.987338
##  [9] 3.989213 3.990059 3.990440 3.990612 3.990690 3.990724 3.990740 3.990747
## [17] 3.990751 3.990752 3.990753 3.990753 3.990753 3.990753 3.990753 3.990753
## [25] 3.990753

sarima(z.ts, 1, 0, 0)

## initial  value 1.394650 
## iter   2 value 1.088119
## iter   3 value 1.088089
## iter   4 value 1.088089
## iter   5 value 1.088089
## iter   5 value 1.088089
## iter   5 value 1.088089
## final  value 1.088089 
## converged
## initial  value 1.087203 
## iter   2 value 1.087196
## iter   3 value 1.087176
## iter   4 value 1.087176
## iter   4 value 1.087176
## iter   4 value 1.087176
## final  value 1.087176 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1    -0.6715 0.0728  -9.2197       0
## xmean  19.8312 0.1776 111.6504       0
## 
## sigma^2 estimated as 8.743891 on 98 degrees of freedom 
##  
## AIC = 5.072228  AICc = 5.073466  BIC = 5.150384 
## 

sarima(z.ts, 2, 0, 0)

## initial  value 1.395908 
## iter   2 value 1.253107
## iter   3 value 1.093980
## iter   4 value 1.079569
## iter   5 value 1.077438
## iter   6 value 1.077335
## iter   7 value 1.077305
## iter   8 value 1.077305
## iter   8 value 1.077305
## final  value 1.077305 
## converged
## initial  value 1.085038 
## iter   2 value 1.084822
## iter   3 value 1.084743
## iter   4 value 1.084718
## iter   5 value 1.084718
## iter   5 value 1.084718
## iter   5 value 1.084718
## final  value 1.084718 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1    -0.6231 0.1002  -6.2215  0.0000
## ar2     0.0704 0.1003   0.7020  0.4843
## xmean  19.8331 0.1906 104.0657  0.0000
## 
## sigma^2 estimated as 8.700424 on 97 degrees of freedom 
##  
## AIC = 5.087314  AICc = 5.089814  BIC = 5.191521 
## 

sarima(z.ts, 1, 0, 1)

## initial  value 1.394650 
## iter   2 value 1.235441
## iter   3 value 1.166536
## iter   4 value 1.125576
## iter   5 value 1.107502
## iter   6 value 1.093580
## iter   7 value 1.088463
## iter   8 value 1.086733
## iter   9 value 1.086515
## iter  10 value 1.086480
## iter  11 value 1.086430
## iter  12 value 1.086419
## iter  13 value 1.086417
## iter  14 value 1.086416
## iter  14 value 1.086416
## iter  14 value 1.086416
## final  value 1.086416 
## converged
## initial  value 1.085073 
## iter   2 value 1.085021
## iter   3 value 1.085011
## iter   4 value 1.085000
## iter   5 value 1.084997
## iter   6 value 1.084996
## iter   7 value 1.084995
## iter   8 value 1.084995
## iter   8 value 1.084995
## final  value 1.084995 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1    -0.7214 0.0976  -7.3942   0.000
## ma1     0.0927 0.1395   0.6644   0.508
## xmean  19.8329 0.1879 105.5291   0.000
## 
## sigma^2 estimated as 8.705361 on 97 degrees of freedom 
##  
## AIC = 5.087867  AICc = 5.090367  BIC = 5.192074 
## 

(예 8-8)

## Example 8.8 단위근 검정
# install.packages("fUnitRoots")
library(astsa) # library for function acf2 & sarima
library(lubridate) # library for function ymd

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(portes) # library for function LjungBox
library(fUnitRoots) # library for function adfTest

z = scan(glue("{dataset_dir}/elecstock.txt"))
stock = ts(z)
ts.plot(stock, ylab = "stock", main = "그림 8-10 주가지수의 시계열그림")

acf2(stock, max.lag = 24, main = "그림 8-10 주가지수의 ACF & PACF")

##      [,1] [,2]  [,3] [,4] [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.96 0.92  0.88 0.85 0.82 0.80  0.78  0.75  0.71  0.67  0.64  0.62  0.59
## PACF 0.96 0.03 -0.07 0.09 0.05 0.03 -0.03 -0.05 -0.12 -0.02  0.07 -0.01  0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.57  0.55  0.53  0.51   0.5  0.48  0.47  0.45  0.43  0.40  0.37
## PACF  0.01  0.01  0.10 -0.07   0.1 -0.01 -0.03 -0.07 -0.03 -0.12 -0.18

LjungBox(stock, lags = seq(6, 24, 6))

##  lags statistic df p-value
##     6  645.8249  6       0
##    12 1074.2322 12       0
##    18 1347.4446 18       0
##    24 1533.4412 24       0

# ADF 검정
adfTest(stock, lags = 0, type = "c") # c : case2 ; 절편이 있고 추세는 없는 경우

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -1.6769
##   P VALUE:
##     0.434 
## 
## Description:
##  Mon Jul  7 03:01:48 2025 by user: PC

adfTest(stock, lags = 1, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -1.5572
##   P VALUE:
##     0.4784 
## 
## Description:
##  Mon Jul  7 03:01:48 2025 by user: PC

adfTest(stock, lags = 2, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -1.629
##   P VALUE:
##     0.4517 
## 
## Description:
##  Mon Jul  7 03:01:48 2025 by user: PC

# function adf.test 를 이용할 수도 있음
# library(tseries) # library for function adf.test & pp.test
# adf.test(stock) # ADF 검정
# pp.test(stock) # PP 검정

dstock = diff(stock, lag = 1)
ts.plot(dstock, ylab = "diff(stock)",
        main = "그림 8-11 차분된 주가지수") ; abline(h = 0)

acf2(dstock, max.lag = 24)

##       [,1] [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.03 0.04 -0.13 -0.08    0 0.05 0.09 0.12 0.00 -0.09  0.00 -0.05  0.01
## PACF -0.03 0.04 -0.13 -0.09    0 0.04 0.07 0.11 0.02 -0.07  0.03 -0.02 -0.03
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.04 -0.11  0.02 -0.09  0.04  0.02  0.06  0.01  0.09  0.12 -0.08
## PACF  0.02 -0.14 -0.01 -0.07  0.02  0.02  0.05  0.03  0.11  0.19 -0.05

LjungBox(dstock, lags = seq(6, 24, 6))

##  lags statistic df   p-value
##     6  4.039155  6 0.6713775
##    12  8.505716 12 0.7444679
##    18 12.237132 18 0.8347560
##    24 17.384384 24 0.8318117

fit = arima(stock, order = c(1, 0, 0), method = "CSS"); fit

## 
## Call:
## arima(x = stock, order = c(1, 0, 0), method = "CSS")
## 
## Coefficients:
##          ar1  intercept
##       0.9603   120.1587
## s.e.  0.0234    13.9993
## 
## sigma^2 estimated as 41.37:  part log likelihood = -442.83

acf2(fit$residuals)

##       [,1] [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.01 0.06 -0.11 -0.07 0.02 0.06 0.10 0.13 0.02 -0.07  0.02 -0.03  0.02
## PACF -0.01 0.06 -0.11 -0.07 0.03 0.06 0.09 0.13 0.03 -0.06  0.05  0.00 -0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## ACF   0.05 -0.10  0.03 -0.07  0.05  0.04  0.07  0.03  0.10
## PACF  0.03 -0.13  0.01 -0.05  0.04  0.04  0.06  0.04  0.12

fit1 = arima(stock, order = c(0,1,0));
summary(fit1)

## 
## Call:
## arima(x = stock, order = c(0, 1, 0))
## 
## 
## sigma^2 estimated as 42.26:  log likelihood = -440.98,  aic = 883.95
## 
## Training set error measures:
##                       ME    RMSE      MAE        MPE     MAPE      MASE
## Training set -0.08043704 6.47675 4.897341 -0.2029409 4.084668 0.9928043
##                     ACF1
## Training set -0.03158845

acf2(fit1$residuals)

##       [,1] [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.03 0.04 -0.13 -0.08    0 0.05 0.09 0.12 0.00 -0.09  0.00 -0.05  0.01
## PACF -0.03 0.04 -0.13 -0.09    0 0.04 0.07 0.11 0.02 -0.07  0.03 -0.02 -0.03
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## ACF   0.04 -0.11  0.02 -0.09  0.04  0.02  0.06  0.01  0.09
## PACF  0.02 -0.14 -0.01 -0.07  0.02  0.02  0.05  0.03  0.11

(예 8-9)

## Example 8.9 : Female Worker
library(astsa) # library for function acf2 & sarima
library(lubridate) # library for function ymd
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest
library(fUnitRoots) # library for function adfTest
library(ggplot2) # library for ggplot

z = scan(glue("{dataset_dir}/female.txt"))
female = ts(z)
ts.plot(female, ylab = "female", main = "그림 8-9 여성 근로자")

date = ymd("821201") + months(1:length(female) - 1)
femaledf = data.frame(female, date)
acf2(female, max.lag = 24, main = "그림 8-9 여성 근로자의 ACF & PACF")

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.97  0.93  0.90  0.87  0.83  0.80 0.76 0.73 0.70  0.68  0.66  0.63  0.61
## PACF 0.97 -0.03 -0.02 -0.01 -0.03 -0.03 0.00 0.03 0.03  0.00  0.04  0.01 -0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.59  0.57  0.55  0.52  0.49  0.46  0.43  0.40  0.37  0.35  0.32
## PACF -0.01 -0.02 -0.02 -0.11 -0.03  0.00 -0.01  0.04 -0.03  0.00 -0.04

LjungBox(female, lags = seq(6, 24, 6))

##  lags statistic df p-value
##     6  532.5488  6       0
##    12  882.5667 12       0
##    18 1119.6281 18       0
##    24 1243.7043 24       0

adfTest(female, lags = 0, type = "ct") # ct : case3 ; 추세가 있는 경우 

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -2.1323
##   P VALUE:
##     0.5217 
## 
## Description:
##  Mon Jul  7 03:01:49 2025 by user: PC

adfTest(female, lags = 1, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.4926
##   P VALUE:
##     0.3723 
## 
## Description:
##  Mon Jul  7 03:01:49 2025 by user: PC

adfTest(female, lags = 2, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.3453
##   P VALUE:
##     0.4334 
## 
## Description:
##  Mon Jul  7 03:01:49 2025 by user: PC

fit1 = lm(female ~ date, data = femaledf) # 선형모형 적합
coefficients(fit1) ; residual = fit1$residuals

##  (Intercept)         date 
## -441.0669259    0.1338777

resdf = data.frame(date, residual)

ggplot(aes(x = date, y = residual), data = resdf) + geom_line()

dfemale = diff(female, lag = 1)
wbar = mean(dfemale) ; gamma_0_hat = var(dfemale)
acf2(dfemale, max.lag = 24) # SACF들은 모두 유의하지 않으므로 S_wbar 계산에 사용하지 않아도 된다.

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF   0.1 -0.07 -0.06 -0.12 -0.02 -0.09 -0.01 -0.14 -0.03  0.02 -0.03  0.07
## PACF  0.1 -0.08 -0.05 -0.11  0.00 -0.11  0.00 -0.18 -0.01 -0.03 -0.06  0.03
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.02  0.01  0.02  0.06 -0.02 -0.06 -0.05 -0.16  0.08  0.00  0.02  0.15
## PACF -0.05 -0.02  0.00  0.05 -0.06 -0.03 -0.07 -0.15  0.08 -0.07  0.02  0.13

S_wbar = sqrt(gamma_0_hat / length(dfemale))
t_0 = wbar / S_wbar ; t_alpha = qt(0.05, (length(dfemale)-2), lower.tail = FALSE)
wbar ; S_wbar ; t_0 ; t_alpha

## [1] 4.214953

## [1] 1.098164

## [1] 3.838181

## [1] 1.659495

female

## Time Series:
## Start = 1 
## End = 108 
## Frequency = 1 
##   [1] 216 223 229 235 227 236 232 234 237 238 244 270 278 277 274 285 262 243
##  [19] 241 238 259 267 272 269 262 276 297 327 341 323 317 322 337 330 329 342
##  [37] 349 346 345 354 345 354 343 355 362 366 367 348 355 362 373 370 375 395
##  [55] 405 408 399 402 403 402 398 408 427 427 433 442 397 416 425 424 434 434
##  [73] 432 440 469 486 495 506 514 529 528 529 521 537 537 544 557 556 561 569
##  [91] 578 589 570 563 566 571 574 569 596 610 627 655 670 663 664 668 664 667

mean(dfemale)

## [1] 4.214953

ts.plot(dfemale, ylab = "diff(female)",
        main = "그림 8-10 차분된 여성 근로자") ; abline(h = 0)

acf2(dfemale, max.lag = 24, main = "그림 8-10 ACF & PACF")

##      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF   0.1 -0.07 -0.06 -0.12 -0.02 -0.09 -0.01 -0.14 -0.03  0.02 -0.03  0.07
## PACF  0.1 -0.08 -0.05 -0.11  0.00 -0.11  0.00 -0.18 -0.01 -0.03 -0.06  0.03
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.02  0.01  0.02  0.06 -0.02 -0.06 -0.05 -0.16  0.08  0.00  0.02  0.15
## PACF -0.05 -0.02  0.00  0.05 -0.06 -0.03 -0.07 -0.15  0.08 -0.07  0.02  0.13

fit2 = arima(female, order = c(0,1,0), method = "ML")
LjungBox(fit2, lags = seq(6, 24, 6))

##  lags statistic df   p-value
##     6  4.455656  6 0.6152616
##    12  7.814954 12 0.7994187
##    18  9.039620 18 0.9588178
##    24 17.333097 24 0.8341120

연습문제

8.1

source(glue(tsa_function_dir, "/ex8_1.R", sep = ""), echo = F)

(a)

sacf_vec = c(-0.48, 0.04, -0.06, 0.14, -0.22, 0.19, -0.10, -0.02, 0.09, 0.03, -0.12, 0.09, 0.03)
spacf_vec = c(-0.48, -0.24, -0.21, 0.01, -0.20, -0.01, -0.05, -0.15, 0.04, 0.06, -0.04, 0.02, 0.06)
n = 100 ; z_bar = 24.75 ; z_sd = 9.27


ex8_1_a = ex8_1(sacf_vec, spacf_vec, n = n, z_bar = z_bar, z_sd = z_sd )

## |spacf[k]| = -0.48, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## |spacf[k]| = -0.24, crit = 0.196
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(2) is acceptable
## |spacf[k]| = -0.21, crit = 0.196
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(3) is acceptable
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(4) is unacceptable
## |spacf[k]| = -0.2, crit = 0.196
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(5) is acceptable
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(6) is unacceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## H_0 : phi_1111 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(11) is unacceptable
## H_0 : phi_1212 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(12) is unacceptable
## H_0 : phi_1313 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(13) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(2) is unacceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(3) is unacceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(4) is unacceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(5) is unacceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(6) is unacceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(7) is unacceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(8) is unacceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(9) is unacceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## H_0 : rho_11 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(11) is unacceptable
## H_0 : rho_12 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(12) is unacceptable
## H_0 : rho_13 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(13) is unacceptable
## MA(1)
## 적절한 모형 : AR(3), MA(1) ; 변수의 절약에 의해 MA(1) 선택
## s_wbar = 0.1854
## t_0 = 133.495145631068, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

(b)

sacf_vec = c(0.59, 0.44, 0.33, 0.23, 0.24, 0.16, 0.05, 0.01, -0.03, -0.11, -0.08, -0.07, 0.01)
spacf_vec = c(0.59, 0.13, 0.04, 0.02, 0.12, -0.05, -0.12, -0.03, -0.01, -0.13, 0.06, 0.04, 0.13)
n = 100 ; z_bar = 25.56 ; z_sd = 3.85


ex8_1_b = ex8_1(sacf_vec, spacf_vec, n = n, z_bar = z_bar, z_sd = z_sd )

## |spacf[k]| = 0.59, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(2) is unacceptable
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(3) is unacceptable
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(4) is unacceptable
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(5) is unacceptable
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(6) is unacceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## H_0 : phi_1111 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(11) is unacceptable
## H_0 : phi_1212 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(12) is unacceptable
## H_0 : phi_1313 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(13) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(2) is acceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(3) is acceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(4) is unacceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(5) is unacceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(6) is unacceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(7) is unacceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(8) is unacceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(9) is unacceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## H_0 : rho_11 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(11) is unacceptable
## H_0 : rho_12 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(12) is unacceptable
## H_0 : rho_13 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(13) is unacceptable
## AR(1)
## 적절한 모형 : AR(1), MA(3) ; 변수의 절약에 의해 AR(1) 선택
## s_wbar = 0.568445687818986
## t_0 = 44.9647179101115, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

(c)

sacf_vec = c(0.97, 0.94, 0.91, 0.88, 0.85, 0.82, 0.79, 0.76, 0.73, 0.70, 0.67, 0.64, 0.61)
spacf_vec = c(0.97, -0.02, -0.01, -0.02, -0.02, -0.01, -0.02, -0.02, -0.02, -0.01, -0.01, -0.02, -0.01)
n = 100 ; z_bar = 548.52 ; z_sd = 149.90


ex8_1_c = ex8_1(sacf_vec, spacf_vec, n = n, z_bar = z_bar, z_sd = z_sd )

## |spacf[k]| = 0.97, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(2) is unacceptable
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(3) is unacceptable
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(4) is unacceptable
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(5) is unacceptable
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(6) is unacceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## H_0 : phi_1111 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(11) is unacceptable
## H_0 : phi_1212 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(12) is unacceptable
## H_0 : phi_1313 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(13) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(2) is acceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(3) is acceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(4) is acceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(5) is acceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(6) is acceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(7) is acceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(8) is acceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(9) is acceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## H_0 : rho_11 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(11) is unacceptable
## H_0 : rho_12 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(12) is unacceptable
## H_0 : rho_13 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(13) is unacceptable
## AR(1)
## 적절한 모형 : AR(1), MA(9) ; 변수의 절약에 의해 AR(1) 선택
## s_wbar = 25.7024958710239
## t_0 = 21.3411181059028, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

(d)

sacf_vec = c(-0.44, -0.05, -0.01, -0.03, 0.12, -0.15, 0.15, -0.04, -0.10, 0.09, 0.08, -0.07, 0.06)
spacf_vec = c(-0.44, -0.31, -0.25, -0.25, -0.07, -0.21, -0.01, 0.02, -0.09, -0.02, 0.03, -0.02, 0.01)
n = 100 ; z_bar = 19.02 ; z_sd = 1.351


ex8_1_d = ex8_1(sacf_vec, spacf_vec, n = n, z_bar = z_bar, z_sd = z_sd, diff = TRUE)

## |spacf[k]| = -0.44, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## |spacf[k]| = -0.31, crit = 0.196
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(2) is acceptable
## |spacf[k]| = -0.25, crit = 0.196
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(3) is acceptable
## |spacf[k]| = -0.25, crit = 0.196
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(4) is acceptable
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(5) is unacceptable
## |spacf[k]| = -0.21, crit = 0.196
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(6) is acceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## H_0 : phi_1111 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(11) is unacceptable
## H_0 : phi_1212 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(12) is unacceptable
## H_0 : phi_1313 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(13) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(2) is unacceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(3) is unacceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(4) is unacceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(5) is unacceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(6) is unacceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(7) is unacceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(8) is unacceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(9) is unacceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## H_0 : rho_11 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(11) is unacceptable
## H_0 : rho_12 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(12) is unacceptable
## H_0 : rho_13 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(13) is unacceptable
## MA(1)
## 적절한 모형 : AR(4), MA(1) ; 변수의 절약에 의해 MA(1) 선택
## t_0 = 14.0784603997039, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

8.2

(a)

library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest

ex8_2a = scan(glue("{dataset_dir}/ex8_2a.txt"))
ex8_2a = ts(ex8_2a) ; ex8_2a = log(ex8_2a)
ts.plot(ex8_2a)

acf2(ex8_2a, max.lag = 24, main = "ex8_2a : SACF, SPACF")

##       [,1]  [,2]  [,3]  [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12]
## ACF  -0.39 -0.19  0.23 -0.20  0.04 0.16 -0.10 -0.08 0.14  0.01 -0.19  0.06
## PACF -0.39 -0.40 -0.05 -0.24 -0.13 0.02  0.02 -0.11 0.02  0.12 -0.11 -0.17
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.05 -0.05  0.14 -0.09 -0.02  0.04 -0.03 -0.01  0.13 -0.10 -0.07  0.14
## PACF -0.09 -0.05  0.07 -0.03  0.07  0.02 -0.03 -0.06  0.16  0.02 -0.12  0.01

fit1 = arima(ex8_2a, order = c(2,0,0)) # 절편 있는 AR(2) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1       -0.546347   0.091487  -5.9719 2.345e-09 ***
## ar2       -0.400105   0.091087  -4.3926 1.120e-05 ***
## intercept  5.822832   0.015971 364.5810 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/100

## [1] 0.09529587

# fit2 = arima(ex8_2a, order = c(2,0,0), include.mean = F) # 절편 없는 AR(2) 모형 적합
# summary(fit2)

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4] [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.02 -0.12 -0.15 -0.15 0.09 0.09 -0.04 -0.03 0.05 -0.04 -0.20 -0.05  0.09
## PACF -0.02 -0.12 -0.15 -0.18 0.04 0.03 -0.06 -0.03 0.08 -0.04 -0.24 -0.08  0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.08  0.16 -0.07 -0.06 -0.05  0.02  0.07   0.1 -0.06 -0.07  0.04
## PACF -0.01  0.09  0.00  0.05 -0.05  0.01  0.06   0.1 -0.10 -0.05  0.10

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df    p-value
##     6  7.783067  4 0.09985584
##    12 13.410133 10 0.20163591
##    18 19.456914 16 0.24568230
##    24 22.672358 22 0.42042344

sarima(ex8_2a, 2,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value -0.993182 
## iter   2 value -1.134944
## iter   3 value -1.144477
## iter   4 value -1.164069
## iter   5 value -1.165861
## iter   6 value -1.166081
## iter   7 value -1.166082
## iter   7 value -1.166082
## iter   7 value -1.166082
## final  value -1.166082 
## converged
## initial  value -1.172776 
## iter   2 value -1.172782
## iter   3 value -1.172799
## iter   4 value -1.172811
## iter   5 value -1.172814
## iter   6 value -1.172814
## iter   6 value -1.172814
## iter   6 value -1.172814
## final  value -1.172814 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1    -0.5463 0.0915  -5.9719       0
## ar2    -0.4001 0.0911  -4.3926       0
## xmean   5.8228 0.0160 364.5810       0
## 
## sigma^2 estimated as 0.09529587 on 97 degrees of freedom 
##  
## AIC = 0.5722493  AICc = 0.5747493  BIC = 0.6764561 
## 

sarima.for(ex8_2a, 6, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 5.915245 5.944070 5.719620 5.830714 5.859822 5.799469
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 0.3087003 0.3517687 0.3531645 0.3631602 0.3647187 0.3650460

sarima.for(ex8_2a, 12, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 5.915245 5.944070 5.719620 5.830714 5.859822 5.799469 5.820797 5.833292
##  [9] 5.817932 5.821325 5.825617 5.821914
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 0.3087003 0.3517687 0.3531645 0.3631602 0.3647187 0.3650460 0.3657035
##  [8] 0.3657494 0.3657923 0.3658319 0.3658326 0.3658368

sarima.for(ex8_2a, 18, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 118 
## Frequency = 1 
##  [1] 5.915245 5.944070 5.719620 5.830714 5.859822 5.799469 5.820797 5.833292
##  [9] 5.817932 5.821325 5.825617 5.821914 5.822220 5.823534 5.822694 5.822627
## [17] 5.823000 5.822823
## 
## $se
## Time Series:
## Start = 101 
## End = 118 
## Frequency = 1 
##  [1] 0.3087003 0.3517687 0.3531645 0.3631602 0.3647187 0.3650460 0.3657035
##  [8] 0.3657494 0.3657923 0.3658319 0.3658326 0.3658368 0.3658390 0.3658390
## [15] 0.3658393 0.3658394 0.3658394 0.3658395

sarima.for(ex8_2a, 24, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 124 
## Frequency = 1 
##  [1] 5.915245 5.944070 5.719620 5.830714 5.859822 5.799469 5.820797 5.833292
##  [9] 5.817932 5.821325 5.825617 5.821914 5.822220 5.823534 5.822694 5.822627
## [17] 5.823000 5.822823 5.822770 5.822870 5.822837 5.822815 5.822840 5.822835
## 
## $se
## Time Series:
## Start = 101 
## End = 124 
## Frequency = 1 
##  [1] 0.3087003 0.3517687 0.3531645 0.3631602 0.3647187 0.3650460 0.3657035
##  [8] 0.3657494 0.3657923 0.3658319 0.3658326 0.3658368 0.3658390 0.3658390
## [15] 0.3658393 0.3658394 0.3658394 0.3658395 0.3658395 0.3658395 0.3658395
## [22] 0.3658395 0.3658395 0.3658395

(b)

library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest

ex8_2b = scan(glue("{dataset_dir}/ex8_2b.txt"))
ex8_2b = ts(ex8_2b) ; n = length(ex8_2b) ; ex8_2b = log(ex8_2b)
ts.plot(ex8_2b) 

acf2(ex8_2b, max.lag = 24, main = "ex8_2b : SACF, SPACF")

##      [,1] [,2]  [,3]  [,4] [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.59 0.46  0.25  0.11 0.05 0.04 -0.03 -0.04 -0.02 -0.05 -0.07 -0.06 -0.06
## PACF 0.59 0.17 -0.13 -0.07 0.03 0.05 -0.08 -0.03  0.07 -0.04 -0.07  0.02  0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.03  0.09  0.07  0.03 -0.05 -0.19 -0.17 -0.21 -0.21 -0.19 -0.19
## PACF  0.01  0.16 -0.05 -0.11 -0.10 -0.18  0.08 -0.06 -0.08  0.00 -0.07

fit1 = arima(ex8_2b, order = c(1,0,0)) # 절편 있는 AR(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1       0.6215577  0.0809138   7.6817  1.57e-14 ***
## intercept 4.6088498  0.0082792 556.6750 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 0.001010076

# fit2 = arima(ex8_2a, order = c(2,0,0), include.mean = F) # 절편 없는 AR(2) 모형 적합
# summary(fit2)

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##      [,1] [,2]  [,3]  [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.1 0.21 -0.03 -0.01 -0.05 0.07 -0.04  0.03 0.05  0.00 -0.04  0.00 -0.05
## PACF -0.1 0.20  0.00 -0.06 -0.05 0.08 -0.01 -0.01 0.06  0.01 -0.07 -0.02 -0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.08  0.11  0.06  0.06  0.06 -0.23 -0.01 -0.08 -0.10  0.01 -0.01
## PACF -0.09  0.10  0.12  0.03  0.02 -0.25 -0.04  0.01 -0.11  0.01  0.01

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  6.496386  5 0.2608674
##    12  7.214314 11 0.7814718
##    18 10.997815 17 0.8566775
##    24 19.834540 23 0.6518714

sarima(ex8_2b, 1,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value -3.212757 
## iter   2 value -3.442672
## iter   3 value -3.443404
## iter   4 value -3.443908
## iter   5 value -3.443947
## iter   6 value -3.443998
## iter   7 value -3.443999
## iter   8 value -3.443999
## iter   8 value -3.443999
## iter   8 value -3.443999
## final  value -3.443999 
## converged
## initial  value -3.446384 
## iter   2 value -3.446413
## iter   3 value -3.446422
## iter   4 value -3.446423
## iter   5 value -3.446423
## iter   5 value -3.446423
## iter   5 value -3.446423
## final  value -3.446423 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1     0.6216 0.0809   7.6817       0
## xmean   4.6088 0.0083 556.6750       0
## 
## sigma^2 estimated as 0.001010076 on 98 degrees of freedom 
##  
## AIC = -3.99497  AICc = -3.993733  BIC = -3.916815 
## 

sarima.for(ex8_2b, 6, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 4.547943 4.570993 4.585319 4.594224 4.599759 4.603199
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 0.03178169 0.03742061 0.03938350 0.04011612 0.04039560 0.04050305

sarima.for(ex8_2b, 12, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 4.547943 4.570993 4.585319 4.594224 4.599759 4.603199 4.605338 4.606667
##  [9] 4.607493 4.608006 4.608326 4.608524
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 0.03178169 0.03742061 0.03938350 0.04011612 0.04039560 0.04050305
##  [7] 0.04054449 0.04056049 0.04056666 0.04056905 0.04056997 0.04057033

(c)

library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest

ex8_2c = scan(glue("{dataset_dir}/ex8_2c.txt"))
ex8_2c = ts(ex8_2c) ; n = length(ex8_2c)
ts.plot(ex8_2c) 

acf2(ex8_2c, max.lag = 24, main = "ex8_2c : SACF, SPACF")

##       [,1]  [,2] [,3]  [,4]  [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.41 -0.15 0.37 -0.16 -0.01 0.21 -0.15 -0.10  0.21 -0.06 -0.21  0.19
## PACF -0.41 -0.39 0.16  0.07  0.09 0.20  0.06 -0.16 -0.04  0.04 -0.18 -0.08
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.09 -0.10  0.13 -0.16 -0.01  0.03 -0.09 -0.09  0.15 -0.09 -0.09  0.10
## PACF -0.09 -0.04 -0.01 -0.11 -0.03 -0.13 -0.14 -0.22  0.04  0.03  0.00 -0.04

adfTest(ex8_2c, lags = 0, type = "c") # c : case2 ; 절편이 있고 추세는 없는 경우

## Warning in adfTest(ex8_2c, lags = 0, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -15.2737
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:52 2025 by user: PC

adfTest(ex8_2c, lags = 1, type = "c")

## Warning in adfTest(ex8_2c, lags = 1, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -12.4898
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:52 2025 by user: PC

adfTest(ex8_2c, lags = 2, type = "c")

## Warning in adfTest(ex8_2c, lags = 2, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -6.3759
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:52 2025 by user: PC

fit1 = arima(ex8_2c, order = c(2,0,0)) # 절편 있는 AR(2) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##             Estimate Std. Error   z value  Pr(>|z|)    
## ar1        -0.574883   0.091825   -6.2606 3.834e-10 ***
## ar2        -0.391975   0.091842   -4.2679 1.973e-05 ***
## intercept 250.004985   0.043637 5729.2594 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 0.7263745

# fit2 = arima(ex8_2c, order = c(0,0,1)) # 절편 있는 MA(1) 모형 적합
# coeftest(fit2)
# sum(fit2$residuals^2)/n

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##      [,1] [,2] [,3] [,4] [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.06 0.09 0.18 0.09 0.15 0.12 -0.08 -0.05 0.06 -0.07 -0.22 -0.02 -0.11
## PACF 0.06 0.08 0.17 0.07 0.12 0.07 -0.14 -0.12 0.03 -0.06 -0.23  0.01 -0.03
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.13 -0.06 -0.24 -0.15 -0.20 -0.17 -0.13  0.03 -0.10 -0.08  0.02
## PACF -0.07 -0.02 -0.13 -0.08 -0.22 -0.11 -0.04  0.13 -0.03 -0.02  0.02

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

# ts.plot(fit2$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)
# acf2(fit2$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")
# 
# qqnorm(resid(fit2), main = "적합된 모형의 잔차의 정규성검정")
# qqline(resid(fit2), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df     p-value
##     6  9.346302  4 0.053003023
##    12 16.992870 10 0.074521899
##    18 35.206429 16 0.003723669
##    24 43.284090 22 0.004353274

sarima(ex8_2c, 2,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 0.024881 
## iter   2 value -0.120118
## iter   3 value -0.132751
## iter   4 value -0.155603
## iter   5 value -0.155712
## iter   6 value -0.155724
## iter   7 value -0.155725
## iter   8 value -0.155725
## iter   8 value -0.155725
## final  value -0.155725 
## converged
## initial  value -0.157204 
## iter   2 value -0.157235
## iter   3 value -0.157239
## iter   4 value -0.157242
## iter   5 value -0.157242
## iter   5 value -0.157242
## iter   5 value -0.157242
## final  value -0.157242 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE   t.value p.value
## ar1    -0.5749 0.0918   -6.2606       0
## ar2    -0.3920 0.0918   -4.2679       0
## xmean 250.0050 0.0436 5729.2594       0
## 
## sigma^2 estimated as 0.7263745 on 97 degrees of freedom 
##  
## AIC = 2.603394  AICc = 2.605894  BIC = 2.707601 
## 

sarima.for(ex8_2c, 6, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 250.1927 250.2910 249.7670 250.0297 250.0841 249.9498
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 0.8522760 0.9830739 0.9844695 1.0092287 1.0149045 1.0152241

sarima.for(ex8_2c, 12, 2, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 250.1927 250.2910 249.7670 250.0297 250.0841 249.9498 250.0057 250.0262
##  [9] 249.9925 250.0038 250.0105 250.0022
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 0.8522760 0.9830739 0.9844695 1.0092287 1.0149045 1.0152241 1.0168037
##  [8] 1.0170545 1.0170963 1.0171948 1.0172049 1.0172092

(d)

library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest

ex8_2d = scan(glue("{dataset_dir}/ex8_2d.txt"))
ex8_2d = ts(ex8_2d) ; n = length(ex8_2d)
ts.plot(ex8_2d) 

acf2(ex8_2d, max.lag = 24, main = "ex8_2d : SACF, SPACF")

##      [,1] [,2]  [,3] [,4] [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.67 0.53  0.36  0.2 0.14 0.10  0.02 -0.04 -0.03 -0.06 -0.07 -0.04 -0.06
## PACF 0.67 0.15 -0.09 -0.1 0.04 0.05 -0.10 -0.09  0.09 -0.02 -0.05  0.03 -0.03
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.03  0.05  0.02 -0.03 -0.12 -0.20 -0.22 -0.24 -0.27 -0.26 -0.28
## PACF  0.06  0.10 -0.11 -0.13 -0.12 -0.08  0.02 -0.08 -0.09  0.00 -0.08

adfTest(ex8_2c, lags = 0, type = "ct") # c : case3 ; 절편이 있고 추세도 있는 경우

## Warning in adfTest(ex8_2c, lags = 0, type = "ct"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -15.2139
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:53 2025 by user: PC

adfTest(ex8_2c, lags = 1, type = "ct")

## Warning in adfTest(ex8_2c, lags = 1, type = "ct"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -12.4421
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:53 2025 by user: PC

adfTest(ex8_2c, lags = 2, type = "ct")

## Warning in adfTest(ex8_2c, lags = 2, type = "ct"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -6.3641
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:53 2025 by user: PC

# fit1 = arima(ex8_2c, order = c(2,0,0)) # 절편 있는 AR(2) 모형 적합
# coeftest(fit1)
# sum(fit1$residuals^2)/n
fit2 = arima(ex8_2d, order = c(1,0,0)) # 절편 있는 AR(1) 모형 적합
coeftest(fit2)

## 
## z test of coefficients:
## 
##             Estimate Std. Error  z value  Pr(>|z|)    
## ar1         0.723930   0.072868   9.9348 < 2.2e-16 ***
## intercept 210.280499   0.394934 532.4447 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 1.23733

# ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)
# acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")
# 
# qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
# qqline(resid(fit1), col = "red")

ts.plot(fit2$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit2$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##      [,1] [,2]  [,3]  [,4]  [,5] [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.1 0.17 -0.01 -0.02 -0.05 0.09 -0.01 -0.01 0.06 -0.01 -0.06  0.06 -0.09
## PACF -0.1 0.16  0.03 -0.05 -0.06 0.10  0.03 -0.04 0.05  0.01 -0.07  0.03 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.04  0.08  0.08  0.06 -0.01 -0.17 -0.06 -0.03 -0.12  0.02 -0.02
## PACF -0.06  0.09  0.12  0.05 -0.06 -0.20 -0.05  0.02 -0.13 -0.01  0.00

qqnorm(resid(fit2), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit2), col = "red")

LjungBox(fit2, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  5.154844  5 0.3972768
##    12  6.345019 11 0.8493806
##    18  9.591094 17 0.9198767
##    24 15.905343 23 0.8592908

sarima(ex8_2d, 1,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 0.450696 
## iter   2 value 0.114639
## iter   3 value 0.112596
## iter   4 value 0.111502
## iter   5 value 0.111443
## iter   6 value 0.111288
## iter   7 value 0.111288
## iter   8 value 0.111288
## iter   8 value 0.111288
## iter   8 value 0.111288
## final  value 0.111288 
## converged
## initial  value 0.110257 
## iter   2 value 0.110215
## iter   3 value 0.110192
## iter   4 value 0.110191
## iter   5 value 0.110191
## iter   6 value 0.110190
## iter   6 value 0.110190
## final  value 0.110190 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE  t.value p.value
## ar1     0.7239 0.0729   9.9348       0
## xmean 210.2805 0.3949 532.4447       0
## 
## sigma^2 estimated as 1.23733 on 98 degrees of freedom 
##  
## AIC = 3.118258  AICc = 3.119495  BIC = 3.196413 
## 

sarima.for(ex8_2d, 6, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 207.1817 208.0372 208.6565 209.1048 209.4294 209.6644
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 1.112353 1.373238 1.491851 1.550394 1.580208 1.595611

sarima.for(ex8_2d, 12, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 207.1817 208.0372 208.6565 209.1048 209.4294 209.6644 209.8345 209.9576
##  [9] 210.0467 210.1113 210.1580 210.1918
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 1.112353 1.373238 1.491851 1.550394 1.580208 1.595611 1.603624 1.607807
##  [9] 1.609995 1.611141 1.611741 1.612055

8.3

(a)

library(glue)
library(astsa) # library for function acf2 & sarima
library(portes) # library for function LjungBox
library(lmtest) # library for function coeftest
library(fUnitRoots)

z = scan(glue(dataset_dir, "/ex7_5a.txt"))
length(z)

## [1] 100

ex7_5a = ts(z)
ts.plot(ex7_5a, ylab = "ex7_5a", main = "모의실험 자료")

acf(ex7_5a, lag.max = 24, main = "ex7_5a의 SACF")

pacf(ex7_5a, lag.max = 24, main = "ex7_5a의 SPACF")

adfTest(ex7_5a, lags = 0, type = "ct") # c : case3 ; 절편이 있고 추세도 있는 경우

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -1.7648
##   P VALUE:
##     0.6741 
## 
## Description:
##  Mon Jul  7 03:01:54 2025 by user: PC

adfTest(ex7_5a, lags = 1, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.2997
##   P VALUE:
##     0.4525 
## 
## Description:
##  Mon Jul  7 03:01:54 2025 by user: PC

adfTest(ex7_5a, lags = 2, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.1676
##   P VALUE:
##     0.5072 
## 
## Description:
##  Mon Jul  7 03:01:54 2025 by user: PC

adfTest(ex7_5a, lags = 3, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 3
##   STATISTIC:
##     Dickey-Fuller: -1.8975
##   P VALUE:
##     0.6191 
## 
## Description:
##  Mon Jul  7 03:01:54 2025 by user: PC

adfTest(ex7_5a, lags = 4, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 4
##   STATISTIC:
##     Dickey-Fuller: -2.1122
##   P VALUE:
##     0.5302 
## 
## Description:
##  Mon Jul  7 03:01:54 2025 by user: PC

fit1 = arima(ex7_5a, order = c(1,1,0)) # 절편 있는 AR(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 0.767230   0.063152  12.149 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 1.23733

diff_a_1 = diff(ex7_5a, lag = 1)
ts.plot(diff_a_1, ylab = "diff_a_1", main = "1차 차분된 모의실험 자료")

acf(diff_a_1, lag.max = 24, main = "1차 차분된 ex7_5a의 SACF")

pacf(diff_a_1, lag.max = 24, main = "1차 차분된 ex7_5a의 SPACF")

# diff_a_2 = diff(diff_a_1, lag = 1)
# ts.plot(diff_a_2, ylab = "diff_a_2", main = "2차 차분된 모의실험 자료")
# acf(diff_a_2, lag.max = 24, main = "2차 차분된 ex7_5a의 SACF")
# pacf(diff_a_2, lag.max = 24, main = "2차 차분된 ex7_5a의 SPACF")

# fit2 = arima(diff_a_1, order = c(1,0,0)) # 절편 있는 AR(1) 모형 적합
# coeftest(fit2)
# sum(fit2$residuals^2)/n

# fit3 = arima(diff_a_1, order = c(1,0,0), include.mean = FALSE) # 절편 없는 AR(1) 모형 적합
# coeftest(fit3)
# sum(fit3$residuals^2)/n

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##      [,1] [,2]  [,3] [,4]  [,5] [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  0.05 0.04 -0.02 0.01 -0.04 -0.2 -0.05 -0.03 -0.10 -0.13  0.04  0.01  0.20
## PACF 0.05 0.04 -0.03 0.02 -0.04 -0.2 -0.02 -0.01 -0.11 -0.12  0.04 -0.04  0.19
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.04  0.00  0.06 -0.14  0.02 -0.07 -0.11  0.02  0.01  0.00  0.02
## PACF  0.02 -0.08  0.03 -0.16  0.02  0.00 -0.13  0.05  0.06 -0.02  0.03

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  5.071382  5 0.4072309
##    12  8.538104 11 0.6644432
##    18 16.394645 17 0.4960505
##    24 18.600503 23 0.7241524

sarima(ex7_5a, 1,1,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 2.884986 
## iter   2 value 2.450885
## iter   3 value 2.450860
## iter   4 value 2.450811
## iter   5 value 2.450750
## iter   6 value 2.450747
## iter   7 value 2.450747
## iter   7 value 2.450747
## iter   7 value 2.450747
## final  value 2.450747 
## converged
## initial  value 2.453793 
## iter   2 value 2.453718
## iter   3 value 2.453531
## iter   4 value 2.453527
## iter   5 value 2.453526
## iter   5 value 2.453526
## iter   5 value 2.453526
## final  value 2.453526 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ar1        0.7578 0.0642 11.8025  0.0000
## constant   3.6999 4.6599  0.7940  0.4291
## 
## sigma^2 estimated as 134.0786 on 97 degrees of freedom 
##  
## AIC = 7.805534  AICc = 7.806797  BIC = 7.884174 
## 

sarima.for(ex7_5a, 6, 1, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 366.3494 364.4793 363.9582 364.4595 365.7356 367.5987
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 11.57923 23.41702 35.74247 48.00212 59.91788 71.35820

sarima.for(ex7_5a, 12, 1, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 366.3494 364.4793 363.9582 364.4595 365.7356 367.5987 369.9068 372.5519
##  [9] 375.4526 378.5468 381.7877 385.1398
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1]  11.57923  23.41702  35.74247  48.00212  59.91788  71.35820  82.27140
##  [8]  92.64990 102.51002 111.88039 120.79519 129.29014

(b)

library(glue)
source(glue(tsa_function_dir, "/ex8_1.R", sep = ""), echo = F)
z = scan(glue(dataset_dir, "/ex7_5b.txt"))
n = length(z)
ex7_5b = ts(z)
ts.plot(ex7_5b, ylab = "ex7_5b", main = "모의실험 자료")

acf(ex7_5b, lag.max = 24, main = "ex7_5b의 SACF")

pacf(ex7_5b, lag.max = 12, main = "ex7_5b의 SPACF")

adfTest(ex7_5b, lags = 0, type = "c") # ct : case3 ; 절편이 있고 추세도 있는 경우

## Warning in adfTest(ex7_5b, lags = 0, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -4.5711
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:55 2025 by user: PC

adfTest(ex7_5b, lags = 1, type = "c")

## Warning in adfTest(ex7_5b, lags = 1, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -4.2143
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:55 2025 by user: PC

adfTest(ex7_5b, lags = 2, type = "c")

## Warning in adfTest(ex7_5b, lags = 2, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -3.687
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:55 2025 by user: PC

# diff_b_1 = diff(ex7_5b, lag = 1)
# ts.plot(diff_b_1, ylab = "diff_b_1", main = "1차 차분된 모의실험 자료")
# acf(diff_b_1, lag.max = 24, main = "1차 차분된 ex7_5b의 SACF")
# pacf(diff_b_1, lag.max = 24, main = "1차 차분된 ex7_5b의 SPACF")

fit1 = arima(ex7_5b, order = c(1,0,0)) # 절편 있는 AR(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1         0.81376    0.06563  12.399 < 2.2e-16 ***
## intercept 156.94935    2.64981  59.230 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 1.23733

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4] [,5]  [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.11 -0.12  0.02 -0.06 0.01  0.00 0.06 -0.03 0.07 -0.10 -0.03 -0.07  0.01
## PACF -0.11 -0.14 -0.01 -0.07 0.00 -0.02 0.06 -0.02 0.09 -0.09 -0.02 -0.12 -0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.02 -0.08  0.07 -0.15 -0.11  0.09 -0.05  0.02  0.11 -0.07 -0.04
## PACF -0.07 -0.09  0.01 -0.16 -0.17  0.02 -0.09  0.00  0.07 -0.06 -0.03

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  3.218451  5 0.6663481
##    12  6.159962 11 0.8624806
##    18 11.780257 17 0.8132569
##    24 15.538792 23 0.8742903

sarima(ex7_5b, 1,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 2.019977 
## iter   2 value 1.580550
## iter   3 value 1.571418
## iter   4 value 1.566079
## iter   5 value 1.565866
## iter   6 value 1.565539
## iter   7 value 1.565538
## iter   8 value 1.565538
## iter   9 value 1.565538
## iter  10 value 1.565538
## iter  10 value 1.565538
## iter  10 value 1.565538
## final  value 1.565538 
## converged
## initial  value 1.645839 
## iter   2 value 1.637141
## iter   3 value 1.633730
## iter   4 value 1.633003
## iter   5 value 1.632714
## iter   6 value 1.632479
## iter   7 value 1.632394
## iter   8 value 1.632252
## iter   9 value 1.632244
## iter  10 value 1.632229
## iter  11 value 1.632221
## iter  12 value 1.632220
## iter  13 value 1.632219
## iter  14 value 1.632219
## iter  15 value 1.632219
## iter  15 value 1.632219
## iter  15 value 1.632219
## final  value 1.632219 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE t.value p.value
## ar1     0.8138 0.0656 12.3991       0
## xmean 156.9493 2.6498 59.2303       0
## 
## sigma^2 estimated as 25.88296 on 98 degrees of freedom 
##  
## AIC = 6.162315  AICc = 6.163552  BIC = 6.24047 
## 

sarima.for(ex7_5b, 6, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 160.2456 159.6317 159.1321 158.7256 158.3948 158.1256
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 5.087530 6.559172 7.373789 7.866943 8.177158 8.376264

sarima.for(ex7_5b, 12, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 160.2456 159.6317 159.1321 158.7256 158.3948 158.1256 157.9065 157.7283
##  [9] 157.5832 157.4652 157.3691 157.2909
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 5.087530 6.559172 7.373789 7.866943 8.177158 8.376264 8.505547 8.590089
##  [9] 8.645617 8.682193 8.706329 8.722275

dl_b = durbin_levinson(ex7_5b, K = 10)

## SACF_1 = 0.722518827798806
## SACF_2 = 0.532281578023179
## SACF_3 = 0.427532290162719
## SACF_4 = 0.323383931010263
## SACF_5 = 0.24562960939634
## SACF_6 = 0.189032660865494
## SACF_7 = 0.147989620322317
## SACF_8 = 0.0777273968860971
## SACF_9 = 0.0338904280641514
## SACF_10 = -0.0697038279939125
## SPACF_1 = 0.722518827798806
## SPACF_2 = 0.0214410854468669
## SPACF_3 = 0.0747323191891475
## SPACF_4 = -0.0326085254834318
## SPACF_5 = 0.00446068098619576
## SPACF_6 = -0.000593629804165394
## SPACF_7 = 0.00776957156653462
## SPACF_8 = -0.0794484550361212
## SPACF_9 = -0.00495072728919832
## SPACF_10 = -0.174012303110664

sacf = dl_b$SACF ; spacf = dl_b$SPACF
sacf

##  [1]  0.72251883  0.53228158  0.42753229  0.32338393  0.24562961  0.18903266
##  [7]  0.14798962  0.07772740  0.03389043 -0.06970383

spacf

##  [1]  0.7225188278  0.0214410854  0.0747323192 -0.0326085255  0.0044606810
##  [6] -0.0005936298  0.0077695716 -0.0794484550 -0.0049507273 -0.1740123031

z_bar = mean(ex7_5b) ; z_sd = sd(ex7_5b)
ex8_1(sacf = sacf, spacf = spacf, n = n, z_bar = z_bar, z_sd = z_sd, diff = FALSE)

## |spacf[k]| = 0.722518827798806, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(2) is unacceptable
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(3) is unacceptable
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(4) is unacceptable
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(5) is unacceptable
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(6) is unacceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(2) is acceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(3) is acceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(4) is acceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(5) is acceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(6) is unacceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(7) is unacceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(8) is unacceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(9) is unacceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## AR(1)
## 적절한 모형 : AR(1), MA(5) ; 변수의 절약에 의해 AR(1) 선택
## s_wbar = 1.25789406340418
## t_0 = 125.567012831389, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

(c)

library(glue)
z = scan(glue(dataset_dir, "/ex7_5c.txt"))
n = length(z) ; n

## [1] 100

ex7_5c = ts(z)
ts.plot(ex7_5c, ylab = "ex7_5c", main = "모의실험 자료")

acf(ex7_5c, lag.max = 24, main = "ex7_5c의 SACF")

pacf(ex7_5c, lag.max = 12, main = "ex7_5c의 SPACF")

adfTest(ex7_5c, lags = 0, type = "c") # c : case2 ; 절편이 있고 추세는 없는 경우

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -2.163
##   P VALUE:
##     0.2545 
## 
## Description:
##  Mon Jul  7 03:01:56 2025 by user: PC

adfTest(ex7_5c, lags = 1, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.0184
##   P VALUE:
##     0.03861 
## 
## Description:
##  Mon Jul  7 03:01:56 2025 by user: PC

adfTest(ex7_5c, lags = 2, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.6091
##   P VALUE:
##     0.09538 
## 
## Description:
##  Mon Jul  7 03:01:56 2025 by user: PC

diff_c_1 = diff(ex7_5c, lag = 1)
ts.plot(diff_c_1, ylab = "diff_c_1", main = "1차 차분된 모의실험 자료")

acf(diff_c_1, lag.max = 24, main = "1차 차분된 ex7_5c의 SACF")

pacf(diff_c_1, lag.max = 24, main = "1차 차분된 ex7_5c의 SPACF")

fit1 = arima(ex7_5c, order = c(0,1,1)) # 차분된 MA(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 0.490465   0.098624  4.9731 6.589e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 101.1674

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.04 -0.07 -0.01 -0.08 0.05 0.02 0.05 0.01 0.01 -0.14 -0.05 -0.09  0.03
## PACF -0.04 -0.07 -0.02 -0.08 0.04 0.01 0.05 0.01 0.03 -0.14 -0.06 -0.12  0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.03 -0.10  0.04 -0.17 -0.13  0.09  0.01  0.03  0.04 -0.11  0.04
## PACF -0.07 -0.11  0.01 -0.17 -0.16  0.05 -0.02 -0.01  0.02 -0.10  0.02

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  1.610159  5 0.9000181
##    12  5.268267 11 0.9174815
##    18 12.559940 17 0.7651106
##    24 15.944937 23 0.8576157

sarima(ex7_5c, 0,1,1) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 2.404331 
## iter   2 value 2.324240
## iter   3 value 2.315202
## iter   4 value 2.313319
## iter   5 value 2.313219
## iter   6 value 2.313218
## iter   6 value 2.313218
## final  value 2.313218 
## converged
## initial  value 2.314064 
## iter   2 value 2.314052
## iter   3 value 2.314052
## iter   4 value 2.314052
## iter   4 value 2.314052
## iter   4 value 2.314052
## final  value 2.314052 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ma1        0.4896 0.0988  4.9562  0.0000
## constant   0.5815 1.5073  0.3858  0.7005
## 
## sigma^2 estimated as 102.037 on 97 degrees of freedom 
##  
## AIC = 7.526587  AICc = 7.527849  BIC = 7.605227 
## 

sarima.for(ex7_5c, 6, 0, 1, 1) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 157.0887 157.6701 158.2516 158.8331 159.4146 159.9960
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 10.10134 18.12317 23.55550 27.95127 31.74407 35.12973

sarima.for(ex7_5c, 12, 0, 1, 1) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 157.0887 157.6701 158.2516 158.8331 159.4146 159.9960 160.5775 161.1590
##  [9] 161.7405 162.3220 162.9034 163.4849
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 10.10134 18.12317 23.55550 27.95127 31.74407 35.12973 38.21662 41.07216
##  [9] 43.74168 46.25740 48.64318 50.91730

dl_c = durbin_levinson(diff_c_1, K = 10)

## SACF_1 = 0.337052005868872
## SACF_2 = -0.0901532330103759
## SACF_3 = -0.0670503016317994
## SACF_4 = -0.0639384797660444
## SACF_5 = 0.026774608744401
## SACF_6 = 0.0575318500113969
## SACF_7 = 0.0567078831710675
## SACF_8 = 0.0421486131051896
## SACF_9 = -0.037331191988209
## SACF_10 = -0.160661576961136
## SPACF_1 = 0.337052005868872
## SPACF_2 = -0.229871637773007
## SPACF_3 = 0.0569343420154281
## SPACF_4 = -0.0923809620347461
## SPACF_5 = 0.0912601835104974
## SPACF_6 = -0.00855560007322475
## SPACF_7 = 0.0598573005602513
## SPACF_8 = 0.00672004405866938
## SPACF_9 = -0.0431155246195279
## SPACF_10 = -0.141952493470491

sacf = dl_c$SACF ; spacf = dl_c$SPACF

z_bar = mean(ex7_5c) ; z_sd = sd(ex7_5c)
ex8_1(sacf = sacf, spacf = spacf, n = n, z_bar = z_bar, z_sd = z_sd, diff = TRUE)

## |spacf[k]| = 0.337052005868872, crit = 0.196
## H_0 : phi_11 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(1) is acceptable
## |spacf[k]| = -0.229871637773007, crit = 0.196
## H_0 : phi_22 = 0 는 유의수준 alpha = 0.05에서 기각된다, AR(2) is acceptable
## H_0 : phi_33 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(3) is unacceptable
## H_0 : phi_44 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(4) is unacceptable
## H_0 : phi_55 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(5) is unacceptable
## H_0 : phi_66 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(6) is unacceptable
## H_0 : phi_77 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(7) is unacceptable
## H_0 : phi_88 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(8) is unacceptable
## H_0 : phi_99 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(9) is unacceptable
## H_0 : phi_1010 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, AR(10) is unacceptable
## 
## H_0 : rho_1 = 0 는 유의수준 alpha = 0.05에서 기각된다, MA(1) is acceptable
## H_0 : rho_2 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(2) is unacceptable
## H_0 : rho_3 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(3) is unacceptable
## H_0 : rho_4 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(4) is unacceptable
## H_0 : rho_5 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(5) is unacceptable
## H_0 : rho_6 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(6) is unacceptable
## H_0 : rho_7 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(7) is unacceptable
## H_0 : rho_8 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(8) is unacceptable
## H_0 : rho_9 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(9) is unacceptable
## H_0 : rho_10 = 0 는 유의수준 alpha = 0.05에서 기각되지 못한다, MA(10) is unacceptable
## MA(1)
## 적절한 모형 : AR(2), MA(1) ; 변수의 절약에 의해 MA(1) 선택
## t_0 = 4.68511731315588, t_alpha = 1.66055121706573
## H_0 : delta = 0 은 유의수준 alpha = 0.05 에서 기각한다, 상수항은 유의하다.

(d)

library(glue)
z = scan(glue(dataset_dir, "/ex7_5d.txt"))
n = length(z)
ex7_5d = ts(z)
ts.plot(ex7_5d, ylab = "ex7_5d", main = "모의실험 자료")

acf(ex7_5d, lag.max = 24, main = "ex7_5d의 SACF")

pacf(ex7_5d, lag.max = 12, main = "ex7_5d의 SPACF")

adfTest(ex7_5d, lags = 0, type = "ct") # ct : case3 ; 절편이 있고 추세도 있는 경우

## Warning in adfTest(ex7_5d, lags = 0, type = "ct"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -4.9988
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:57 2025 by user: PC

adfTest(ex7_5d, lags = 1, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.6051
##   P VALUE:
##     0.03626 
## 
## Description:
##  Mon Jul  7 03:01:57 2025 by user: PC

adfTest(ex7_5d, lags = 2, type = "ct")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.521
##   P VALUE:
##     0.3608 
## 
## Description:
##  Mon Jul  7 03:01:57 2025 by user: PC

diff_d_1 = diff(ex7_5d, lag = 1)
ts.plot(diff_d_1, ylab = "diff_d_1", main = "1차 차분된 모의실험 자료")

acf(diff_d_1, lag.max = 24, main = "1차 차분된 ex7_5d의 SACF")

pacf(diff_d_1, lag.max = 24, main = "1차 차분된 ex7_5d의 SPACF")

fit1 = arima(ex7_5d, order = c(2,1,0)) # 차분된 AR(2) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.356938   0.095784 -3.7265 0.0001942 ***
## ar2 -0.339485   0.096399 -3.5217 0.0004289 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 68.7794

fit2 = arima(ex7_5d, order = c(0,1,2)) # 차분된 MA(2) 모형 적합
coeftest(fit2)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)   
## ma1 -0.33631    0.12716 -2.6448 0.008173 **
## ma2 -0.10829    0.11377 -0.9519 0.341153   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 71.82533

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4]  [,5] [,6]  [,7] [,8]  [,9] [,10] [,11] [,12]
## ACF  -0.13 -0.15 -0.14 -0.11  0.11 0.13 -0.02 0.01 -0.01  0.00 -0.01  0.04
## PACF -0.13 -0.17 -0.20 -0.21 -0.02 0.07 -0.01 0.05  0.06  0.05  0.01  0.06
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.03 -0.03  0.06 -0.05 -0.07  0.02  0.08  0.01  0.09 -0.07 -0.03 -0.03
## PACF  0.00 -0.03  0.05 -0.04 -0.11 -0.05  0.05 -0.02  0.10  0.02  0.05 -0.02

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df    p-value
##     6  10.50492  4 0.03272933
##    12  10.81881 10 0.37180756
##    18  12.41643 16 0.71485057
##    24  15.17529 22 0.85460608

sarima(ex7_5d, 2,1,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 2.172286 
## iter   2 value 2.046251
## iter   3 value 2.037955
## iter   4 value 2.036654
## iter   5 value 2.035613
## iter   6 value 2.034870
## iter   7 value 2.034867
## iter   8 value 2.034864
## iter   8 value 2.034864
## final  value 2.034864 
## converged
## initial  value 2.057985 
## iter   2 value 2.057929
## iter   3 value 2.057529
## iter   4 value 2.057528
## iter   5 value 2.057528
## iter   6 value 2.057528
## iter   6 value 2.057528
## iter   6 value 2.057528
## final  value 2.057528 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ar1       -0.4323 0.0920 -4.6979   0e+00
## ar2       -0.4209 0.0928 -4.5353   0e+00
## constant   1.6442 0.4264  3.8555   2e-04
## 
## sigma^2 estimated as 60.955 on 96 degrees of freedom 
##  
## AIC = 7.033741  AICc = 7.036292  BIC = 7.138594 
## 

sarima.for(ex7_5d, 6, 2, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 183.3406 184.9662 186.7461 188.3394 189.9484 191.6292
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1]  7.807369  8.977767  9.348219 10.515690 11.489462 12.095163

sarima.for(ex7_5d, 12, 2, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 183.3406 184.9662 186.7461 188.3394 189.9484 191.6292 193.2723 194.9015
##  [9] 196.5525 198.2000 199.8398 201.4845
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1]  7.807369  8.977767  9.348219 10.515690 11.489462 12.095163 12.811963
##  [8] 13.543252 14.159353 14.761713 15.365321 15.930937

(e)

library(glue)
z = scan(glue(dataset_dir, "/ex7_5e.txt"))
length(z)

## [1] 100

ex7_5e = ts(z)
ts.plot(ex7_5e, ylab = "ex7_5e", main = "모의실험 자료")

acf(ex7_5e, lag.max = 24, main = "ex7_5e의 SACF")

pacf(ex7_5e, lag.max = 12, main = "ex7_5e의 SPACF")

adfTest(ex7_5e, lags = 0, type = "c") # c : case2 ; 절편이 있고 추세도 있는 경우

## Warning in adfTest(ex7_5e, lags = 0, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -4.0282
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:58 2025 by user: PC

adfTest(ex7_5e, lags = 1, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -3.3024
##   P VALUE:
##     0.01921 
## 
## Description:
##  Mon Jul  7 03:01:58 2025 by user: PC

adfTest(ex7_5e, lags = 2, type = "c")

## Warning in adfTest(ex7_5e, lags = 2, type = "c"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -3.7029
##   P VALUE:
##     0.01 
## 
## Description:
##  Mon Jul  7 03:01:58 2025 by user: PC

fit1 = arima(ex7_5d, order = c(1,0,0)) # AR(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)    
## ar1         0.991321   0.010384 95.4699  < 2e-16 ***
## intercept 103.089114  58.520544  1.7616  0.07814 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 84.3522

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12]
## ACF  -0.33 -0.25  0.13 -0.08  0.03  0.11 -0.08 -0.01 0.03 -0.03 -0.03  0.08
## PACF -0.33 -0.40 -0.16 -0.25 -0.15 -0.01 -0.04  0.00 0.02  0.01 -0.05  0.03
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.05 -0.04  0.09 -0.03 -0.11  0.04  0.07 -0.06  0.06 -0.05  0.01 -0.03
## PACF -0.03 -0.06  0.02  0.01 -0.11 -0.12 -0.03 -0.07  0.02 -0.02  0.07 -0.04

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df      p-value
##     6  21.32740  5 0.0007024203
##    12  22.97128 11 0.0178406612
##    18  26.17639 17 0.0713062735
##    24  28.17450 23 0.2092546261

sarima(ex7_5e, 1,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 0.398876 
## iter   2 value -0.040831
## iter   3 value -0.042521
## iter   4 value -0.043305
## iter   5 value -0.043326
## iter   6 value -0.043371
## iter   7 value -0.043371
## iter   8 value -0.043371
## iter   8 value -0.043371
## iter   8 value -0.043371
## final  value -0.043371 
## converged
## initial  value -0.001073 
## iter   2 value -0.004049
## iter   3 value -0.005960
## iter   4 value -0.006351
## iter   5 value -0.006623
## iter   6 value -0.006809
## iter   7 value -0.006851
## iter   8 value -0.006892
## iter   9 value -0.006892
## iter  10 value -0.006892
## iter  11 value -0.006893
## iter  12 value -0.006893
## iter  13 value -0.006893
## iter  14 value -0.006893
## iter  15 value -0.006893
## iter  16 value -0.006893
## iter  16 value -0.006893
## final  value -0.006893 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE t.value p.value
## ar1     0.7967 0.0655 12.1657       0
## xmean  12.9581 0.4762 27.2128       0
## 
## sigma^2 estimated as 0.9764233 on 98 degrees of freedom 
##  
## AIC = 2.884091  AICc = 2.885328  BIC = 2.962246 
## 

sarima.for(ex7_5e, 6, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 11.39799 11.71509 11.96773 12.16903 12.32940 12.45719
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 0.9881413 1.2634315 1.4105788 1.4964995 1.5485704 1.5807354

sarima.for(ex7_5e, 12, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 11.39799 11.71509 11.96773 12.16903 12.32940 12.45719 12.55900 12.64011
##  [9] 12.70474 12.75623 12.79726 12.82995
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 0.9881413 1.2634315 1.4105788 1.4964995 1.5485704 1.5807354 1.6008185
##  [8] 1.6134376 1.6213972 1.6264298 1.6296165 1.6316362

fit2 = arima(ex7_5d, order = c(1,0,0), include.mean = FALSE, method = "ML") # AR(1) 모형 적합
coeftest(fit2)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 0.9974139  0.0032955  302.66 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 84.43581

ts.plot(fit2$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit2$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12]
## ACF  -0.31 -0.28  0.12 -0.08  0.04  0.11 -0.08 -0.01 0.03 -0.02 -0.02  0.07
## PACF -0.31 -0.42 -0.17 -0.29 -0.16 -0.04 -0.04  0.01 0.03  0.05  0.00  0.09
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF  -0.05 -0.05  0.11 -0.03 -0.09  0.04  0.05 -0.04  0.07 -0.06  0.01 -0.04
## PACF  0.02 -0.04  0.05  0.02 -0.07 -0.08 -0.01 -0.04  0.07  0.02  0.12  0.00

qqnorm(resid(fit2), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit2), col = "red")

LjungBox(fit2, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df    p-value
##     6  21.90363  5 0.00054615
##    12  23.19846 11 0.01656986
##    18  26.57494 17 0.06460340
##    24  28.45057 23 0.19922163

sarima(ex7_5e, 1,0,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value 0.398876 
## iter   2 value -0.040831
## iter   3 value -0.042521
## iter   4 value -0.043305
## iter   5 value -0.043326
## iter   6 value -0.043371
## iter   7 value -0.043371
## iter   8 value -0.043371
## iter   8 value -0.043371
## iter   8 value -0.043371
## final  value -0.043371 
## converged
## initial  value -0.001073 
## iter   2 value -0.004049
## iter   3 value -0.005960
## iter   4 value -0.006351
## iter   5 value -0.006623
## iter   6 value -0.006809
## iter   7 value -0.006851
## iter   8 value -0.006892
## iter   9 value -0.006892
## iter  10 value -0.006892
## iter  11 value -0.006893
## iter  12 value -0.006893
## iter  13 value -0.006893
## iter  14 value -0.006893
## iter  15 value -0.006893
## iter  16 value -0.006893
## iter  16 value -0.006893
## final  value -0.006893 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE t.value p.value
## ar1     0.7967 0.0655 12.1657       0
## xmean  12.9581 0.4762 27.2128       0
## 
## sigma^2 estimated as 0.9764233 on 98 degrees of freedom 
##  
## AIC = 2.884091  AICc = 2.885328  BIC = 2.962246 
## 

sarima.for(ex7_5e, 6, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 11.39799 11.71509 11.96773 12.16903 12.32940 12.45719
## 
## $se
## Time Series:
## Start = 101 
## End = 106 
## Frequency = 1 
## [1] 0.9881413 1.2634315 1.4105788 1.4964995 1.5485704 1.5807354

sarima.for(ex7_5e, 12, 1, 0, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 11.39799 11.71509 11.96773 12.16903 12.32940 12.45719 12.55900 12.64011
##  [9] 12.70474 12.75623 12.79726 12.82995
## 
## $se
## Time Series:
## Start = 101 
## End = 112 
## Frequency = 1 
##  [1] 0.9881413 1.2634315 1.4105788 1.4964995 1.5485704 1.5807354 1.6008185
##  [8] 1.6134376 1.6213972 1.6264298 1.6296165 1.6316362

8.4

library(glue)
z = scan(glue(dataset_dir, "/interest.txt"))
n = length(z)
interest = ts(z)
ts.plot(interest, ylab = "interest", main = "이자율 자료")

acf(interest, lag.max = 24, main = "이자율의 SACF")

pacf(interest, lag.max = 24, main = "이자율의 SPACF")

adfTest(interest, lags = 0, type = "c") # c : case2 ; 절편이 있고 추세도 있는 경우

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 0
##   STATISTIC:
##     Dickey-Fuller: -2.1897
##   P VALUE:
##     0.2432 
## 
## Description:
##  Mon Jul  7 03:01:59 2025 by user: PC

adfTest(interest, lags = 1, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.2586
##   P VALUE:
##     0.2177 
## 
## Description:
##  Mon Jul  7 03:01:59 2025 by user: PC

adfTest(interest, lags = 2, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 2
##   STATISTIC:
##     Dickey-Fuller: -2.3139
##   P VALUE:
##     0.1971 
## 
## Description:
##  Mon Jul  7 03:01:59 2025 by user: PC

diff_int = diff(interest, lag = 1)
ts.plot(diff_int, ylab = "diff_int", main = "1차 차분된 이자율 자료")

acf(diff_int, lag.max = 24, main = "1차 차분된 이자율의 SACF")

pacf(diff_int, lag.max = 24, main = "1차 차분된 이자율의 SPACF")

fit1 = arima(interest, order = c(1,1,0)) # 1차 차분된 AR(1) 모형 적합
coeftest(fit1)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 0.316737   0.075887  4.1738 2.996e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit1$residuals^2)/n

## [1] 0.2553032

fit2 = arima(interest, order = c(0,1,1)) # 1차 차분된 MA(1) 모형 적합
coeftest(fit2)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 0.295422   0.075824  3.8961 9.773e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sum(fit2$residuals^2)/n

## [1] 0.2580358

ts.plot(fit1$residuals, main = "적합된 모형의 잔차 시계열그림"); abline(h = 0)

acf2(fit1$residuals, max.lag = 24, main = "적합된 모형의 잔차 SACF와 SPACF")

##       [,1] [,2] [,3]  [,4]  [,5]  [,6]  [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  -0.03 0.02 0.05 -0.14  0.00 -0.01 -0.11 0.10 0.07  -0.1  0.02  0.12 -0.04
## PACF -0.03 0.02 0.05 -0.13 -0.01 -0.01 -0.10 0.08 0.08  -0.1 -0.02  0.15 -0.01
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## ACF   0.00  0.12 -0.15  0.02 -0.09 -0.07  0.10 -0.02 -0.07  0.03 -0.04
## PACF -0.05  0.15 -0.12 -0.04 -0.06  0.01  0.03 -0.02 -0.04 -0.04 -0.03

qqnorm(resid(fit1), main = "적합된 모형의 잔차의 정규성검정")
qqline(resid(fit1), col = "red")

LjungBox(fit1, lags = seq(6, 24, 6)) # 잔차의 포트맨토 검정

##  lags statistic df   p-value
##     6  3.806045  5 0.5776650
##    12 12.932207 11 0.2977795
##    18 21.755995 17 0.1942859
##    24 26.138732 23 0.2943607

sarima(interest, 1,1,0) # sarima를 이용한 모형 적합을 이용한 경우

## initial  value -0.659850 
## iter   2 value -0.708787
## iter   3 value -0.709020
## iter   4 value -0.709037
## iter   5 value -0.709049
## iter   5 value -0.709049
## iter   5 value -0.709049
## final  value -0.709049 
## converged
## initial  value -0.680849 
## iter   2 value -0.681173
## iter   3 value -0.681294
## iter   4 value -0.681322
## iter   5 value -0.681325
## iter   6 value -0.681325
## iter   6 value -0.681325
## final  value -0.681325 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##          Estimate     SE t.value p.value
## ar1        0.3131 0.0759  4.1260  0.0001
## constant  -0.0465 0.0570 -0.8152  0.4162
## 
## sigma^2 estimated as 0.2558226 on 164 degrees of freedom 
##  
## AIC = 1.511372  AICc = 1.511816  BIC = 1.567613 
## 

sarima.for(interest, 6, 1, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 168 
## End = 173 
## Frequency = 1 
## [1] 11.80738 11.75271 11.70367 11.65639 11.60966 11.56311
## 
## $se
## Time Series:
## Start = 168 
## End = 173 
## Frequency = 1 
## [1] 0.5057891 0.8348003 1.0982983 1.3183393 1.5089356 1.6786869

sarima.for(interest, 12, 1, 1, 0) # 예측값의 시계열 그림

## $pred
## Time Series:
## Start = 168 
## End = 179 
## Frequency = 1 
##  [1] 11.80738 11.75271 11.70367 11.65639 11.60966 11.56311 11.51660 11.47012
##  [9] 11.42364 11.37716 11.33068 11.28420
## 
## $se
## Time Series:
## Start = 168 
## End = 179 
## Frequency = 1 
##  [1] 0.5057891 0.8348003 1.0982983 1.3183393 1.5089356 1.6786869 1.8329740
##  [8] 1.9753018 2.1080582 2.2329404 2.3512004 2.4637909